clear all;
clc;

%% ----------------------- Clinical Data ----------------------------------

p_per_lam=5; % percent model DC fibers / lamina
peract_sens_sub=p_per_lam*[2,2,2,2,2];
peract_sens_epi=p_per_lam*[2,2,3,2,2];
peract_disc_sub=p_per_lam*[8,2,7,3,4];
peract_disc_epi=p_per_lam*[6,4,6,6,5];

nodiam=5;
nopat=5;
diamval=[3,6,9,12,15];
xlab_diam=repmat(diamval,[1,nopat]);

dw=0.65;
temp_base=(-dw/2):dw/(nodiam-1):(dw/2);
w_base=repmat(temp_base,[1,nopat]);
w_base=w_base(:);
temp_trans=repmat(1:nodiam,[nopat,1]);
w_trans=temp_trans(:);
xval_diam=w_base+w_trans;

%% ----------------------- Import Data ------------------------------------

mrg_or_sw=1; % 1 = MRG, 2 = SW

% data format = 5x5
%  - rows: diameters of 3, 6, 9, 12, and 15 um
%  - columns: patients 1-5

if(mrg_or_sw==1)
    cd('mrg_axons');

    perdc_disc_epi=load('perdc_disc_epi.txt');
    perdc_disc_epi=perdc_disc_epi(:);
    perdc_disc_sub=load('perdc_disc_sub.txt');
    perdc_disc_sub=perdc_disc_sub(:);

    perdc_sens_epi=load('perdc_sens_epi.txt');
    perdc_sens_epi=perdc_sens_epi(:);
    perdc_sens_sub=load('perdc_sens_sub.txt');
    perdc_sens_sub=perdc_sens_sub(:);

    perdr_disc_epi=load('perdr_disc_epi.txt');
    perdr_disc_epi=perdr_disc_epi(:);
    perdr_disc_sub=load('perdr_disc_sub.txt');
    perdr_disc_sub=perdr_disc_sub(:);

    perdr_sens_epi=load('perdr_sens_epi.txt');
    perdr_sens_epi=perdr_sens_epi(:);
    perdr_sens_sub=load('perdr_sens_sub.txt');
    perdr_sens_sub=perdr_sens_sub(:);

    cd ..;
else
    cd('sw_axons');

    perdc_disc_epi=load('perdc_disc_epi.txt');
    perdc_disc_epi=perdc_disc_epi(:);
    perdc_disc_sub=load('perdc_disc_sub.txt');
    perdc_disc_sub=perdc_disc_sub(:);

    perdc_sens_epi=load('perdc_sens_epi.txt');
    perdc_sens_epi=perdc_sens_epi(:);
    perdc_sens_sub=load('perdc_sens_sub.txt');
    perdc_sens_sub=perdc_sens_sub(:);

    cd ..;    
end


%% ----------------------- Statistics -------------------------------------

indx_3=1:nodiam:(nodiam*nopat);
indx_6=2:nodiam:(nodiam*nopat);
indx_9=3:nodiam:(nodiam*nopat);
indx_12=4:nodiam:(nodiam*nopat);
indx_15=5:nodiam:(nodiam*nopat);

if(mrg_or_sw==1)
    p1_3=[perdc_sens_sub(indx_3);perdc_sens_epi(indx_3)];
    p2_3=[perdc_disc_sub(indx_3);perdc_disc_epi(indx_3)];
    [h3_dcsens_vs_dcdisc,pval_1]=kstest2(p1_3,p2_3)
    p1_6=[perdc_sens_sub(indx_6);perdc_sens_epi(indx_6)];
    p2_6=[perdc_disc_sub(indx_6);perdc_disc_epi(indx_6)];
    h6_dcsens_vs_dcdisc=kstest2(p1_6,p2_6)
    p1_9=[perdc_sens_sub(indx_9);perdc_sens_epi(indx_9)];
    p2_9=[perdc_disc_sub(indx_9);perdc_disc_epi(indx_9)];
    h9_dcsens_vs_dcdisc=kstest2(p1_9,p2_9)
    p1_12=[perdc_sens_sub(indx_12);perdc_sens_epi(indx_12)];
    p2_12=[perdc_disc_sub(indx_12);perdc_disc_epi(indx_12)];
    h12_dcsens_vs_dcdisc=kstest2(p1_12,p2_12)
    p1_15=[perdc_sens_sub(indx_15);perdc_sens_epi(indx_15)];
    p2_15=[perdc_disc_sub(indx_15);perdc_disc_epi(indx_15)];
    h15_dcsens_vs_dcdisc=kstest2(p1_15,p2_15)

    pp1_3=[perdc_sens_sub(indx_3);perdc_sens_epi(indx_3)];
    pp2_3=[perdr_disc_sub(indx_3);perdr_disc_epi(indx_3)];
    h3_dcsens_vs_drdisc=kstest2(pp1_3,pp2_3)
    pp1_6=[perdc_sens_sub(indx_6);perdc_sens_epi(indx_6)];
    pp2_6=[perdr_disc_sub(indx_6);perdr_disc_epi(indx_6)];
    h6_dcsens_vs_drdisc=kstest2(pp1_6,pp2_6)
    pp1_9=[perdc_sens_sub(indx_9);perdc_sens_epi(indx_9)];
    pp2_9=[perdr_disc_sub(indx_9);perdr_disc_epi(indx_9)];
    h9_dcsens_vs_drdisc=kstest2(pp1_9,pp2_9)
    pp1_12=[perdc_sens_sub(indx_12);perdc_sens_epi(indx_12)];
    pp2_12=[perdr_disc_sub(indx_12);perdr_disc_epi(indx_12)];
    h12_dcsens_vs_drdisc=kstest2(pp1_12,pp2_12)
    pp1_15=[perdc_sens_sub(indx_15);perdc_sens_epi(indx_15)];
    pp2_15=[perdr_disc_sub(indx_15);perdc_disc_epi(indx_15)];
    h15_dcsens_vs_drdisc=kstest2(pp1_15,pp2_15)
end


%% ----------------------- Plotting ---------------------------------------

ms_1=15;
ms_2=12;
lw=5;
fs=28;

dx=(1-dw)/2;
x1=xval_diam(1)-dx;
x2=xval_diam(end)+dx;
y1=-5;
y2=105;

del_xa=(1-dw)/2-0.03;

subplot(2,1,1);
hold on;
plot(xval_diam,perdc_sens_sub,'ks','markersize',ms_1,'linewidth',3);
plot(xval_diam,perdc_sens_epi,'ko','markersize',ms_2,...
    'markerfacecolor','w','linewidth',3);
minix_sub=zeros(nopat,1);
minix_epi=zeros(nopat,1);
for ii=1:nopat
    
    ia=1 + nodiam*(ii-1);
    ib=nodiam*ii;
    [junk,tmp_sub]=min( abs(perdc_sens_sub(ia:ib)-peract_sens_sub(ii)) );
    [junk,tmp_epi]=min( abs(perdc_sens_epi(ia:ib)-peract_sens_epi(ii)) );
    minix_sub(ii)=tmp_sub+nodiam*(ii-1);
    minix_epi(ii)=tmp_epi+nodiam*(ii-1);
    
    xa_sub=ii+[-dw/2-del_xa,-dw/2];
    ya_sub=[peract_sens_sub(ii),peract_sens_sub(ii)];
    
    xa_epi=ii+[dw/2,dw/2+del_xa];
    ya_epi=[peract_sens_epi(ii),peract_sens_epi(ii)];
    
    plot(xa_sub,ya_sub,'k--','linewidth',lw);
    plot(xa_epi,ya_epi,'k-','linewidth',lw);
    plot((ii+0.5)*[1,1],[y1,y2],'k-','color',0.5*[1,1,1],...
        'linewidth',3);
    
end
plot(xval_diam(minix_sub),perdc_sens_sub(minix_sub),...
    'ks','markersize',10,'markerfacecolor','k','linewidth',3);
plot(xval_diam(minix_epi),perdc_sens_epi(minix_epi),...
    'ko','markersize',8,'markerfacecolor','k','linewidth',3); 
hold off;
ylabel({'% DC';'Sensory'},'fontsize',30);
set(gca,'xtick',[],'ytick',0:25:100,'fontsize',fs);
xlim([x1,x2]);
ylim([y1,y2]);

subplot(2,1,2);
hold on;
plot(xval_diam,perdc_disc_sub,'ks','markersize',ms_1,'linewidth',3);
plot(xval_diam,perdc_disc_epi,'ko','markersize',ms_2,...
    'markerfacecolor','w','linewidth',3);
for ii=1:nopat

    ia=1 + nodiam*(ii-1);
    ib=nodiam*ii;
    [junk,tmp_sub]=min( abs(perdc_disc_sub(ia:ib)-peract_disc_sub(ii)) );
    [junk,tmp_epi]=min( abs(perdc_disc_epi(ia:ib)-peract_disc_epi(ii)) );
    minix_sub(ii)=tmp_sub+nodiam*(ii-1);
    minix_epi(ii)=tmp_epi+nodiam*(ii-1);
    
    xa_sub=ii+[-dw/2-del_xa,-dw/2];
    ya_sub=[peract_disc_sub(ii),peract_disc_sub(ii)];
    
    xa_epi=ii+[dw/2,dw/2+del_xa];
    ya_epi=[peract_disc_epi(ii),peract_disc_epi(ii)];
    
    plot(xa_sub,ya_sub,'k--','linewidth',lw);
    plot(xa_epi,ya_epi,'k-','linewidth',lw);
    plot((ii+0.5)*[1,1],[y1,y2],'k-','color',0.5*[1,1,1],...
        'linewidth',3);
    
end
plot(xval_diam(minix_sub),perdc_disc_sub(minix_sub),...
    'ks','markersize',10,'markerfacecolor','k','linewidth',3);
plot(xval_diam(minix_epi),perdc_disc_epi(minix_epi),...
    'ko','markersize',8,'markerfacecolor','k','linewidth',3); 
hold off;
ylabel({'% DC';'Discomfort'},'fontsize',30);
set(gca,'xtick',xval_diam,'xticklabel',xlab_diam,'fontsize',fs,...
    'ytick',0:25:100);
xlim([x1,x2]);
ylim([y1,y2]);

% subplot(3,1,3);
% hold on;
% plot(xval_diam,perdr_sens_sub,'ks','markersize',10,'linewidth',3);
% plot(xval_diam,perdr_sens_epi,'ko','markersize',8,...
%     'markerfacecolor','w','linewidth',3);
% for ii=1:nopat
%     
%     ia=1 + nodiam*(ii-1);
%     ib=nodiam*ii;
%     [junk,tmp_sub]=min( abs(perdr_sens_sub(ia:ib)-peract_sens_sub(ii)) );
%     [junk,tmp_epi]=min( abs(perdr_sens_epi(ia:ib)-peract_sens_epi(ii)) );
%     minix_sub(ii)=tmp_sub+nodiam*(ii-1);
%     minix_epi(ii)=tmp_epi+nodiam*(ii-1);    
%     
%     xa_sub=ii+[-dw/2-del_xa,-dw/2];
%     ya_sub=[peract_sens_sub(ii),peract_sens_sub(ii)];
%     
%     xa_epi=ii+[dw/2,dw/2+del_xa];
%     ya_epi=[peract_sens_epi(ii),peract_sens_epi(ii)];
%     
%     plot(xa_sub,ya_sub,'k--','linewidth',4);
%     plot(xa_epi,ya_epi,'k-','linewidth',4);
%     plot((ii+0.5)*[1,1],[y1,y2],'k-','color',0.4*[1,1,1],...
%         'linewidth',3);
%     
% end
% plot(xval_diam(minix_sub),perdr_disc_sub(minix_sub),...
%     'ks','markersize',10,'markerfacecolor','k','linewidth',3);
% plot(xval_diam(minix_epi),perdr_disc_epi(minix_epi),...
%     'ko','markersize',8,'markerfacecolor','k','linewidth',3); 
% hold off;
% ylabel({'% DR';'Discomfort'},'fontsize',24);
% set(gca,'xtick',xval_diam,'xticklabel',xlab_diam,'fontsize',20);
% xlim([x1,x2]);
% ylim([y1,y2]);

%% ----------------------- Coefficient of Variation -----------------------
% sse = sum squared error

% noaxons=200;
% del_per=1/noaxons*100;
% per=1:del_per:100;
% a_w=1/sum(1./(101-per));

